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METHOD AND APPARATUS FOR PRODUCING AN IMAGE CONTAINING 

DEPTH INFORMATION 

Field of the Invention 
5 This invention relates to a method and apparatus for 
producing an image containing depth information. The 
invention has particular application to thermal images in 
the 8 to 12 micron band of the electromagnetic spectrum. 
However, the method and apparatus is also useful with 
10 images taken in other parts of the electromagnetic 

spectrum (and in particular visible light) , as well as 
electron beam images, and other forms of particle 
radiation images as well as acoustic images. 

15 Background Art 

Thermal images are used to produce a visual image of a 
scene in darkness . In such images , articles in the scene 
which are above ambient temperature can be easily 
visualised. However, one major difficulty with 

20 conventional thermal images is that the images contain 
little or no depth information, so that it is very 
difficult to ascertain whether one part of the image is in 
front of or behind another part of the image. Thus, when 
a thermal image is taken of a scene which, for example, 

2 5 may contain people, the thermal image will clearly show 
each of the individual people, but will not provide 
sufficient information to ascertain the relative positions 
of the people in space. 

30 In many applications, it is necessary that the image does 
include depth information so that the image provides 
sufficient information of the scene to suit a particular 
purpose . 

35 Sxmimary of the Invention 

The object of the invention is to provide a method and 
apparatus of producing images containing improved depth 
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informa-kion so -that a determination can be made as to the 
distance of various parts of the image relative to one 
another and from the location where the image is taken. 



5 The invention may be said to reside in a method of 
producing an image containing depth information, 
comprising: 

detecting a radiation wavefield emanating from 
the scene and taking at least two images of the scene at 
10 different planes relative to the scene to produce a first 
image comprised of a first set of intensity data values 
and a second image comprised of a second set of intensity 
data values; 

determining an intensity variation of data values 
15 in the first set of data values relative to other data 

values in the first set of data values to produce a first 
set of intensity variances, and determining an intensity 
variation of data values in the second set of values 
relative to other data values in the second set of values 
20 to produce a second set of intensity variance data; 

processing the first and second sets of intensity 
variance data to obtain image data of the scene containing 
depth information; and 

coding the image data which have the same depth 
25 information with a code reference to identify the 
different depth information in the image data. 



Thus, by producing the variance data and then processing 
the variance data, the variance data is able to result in 
30 the production of depth information in the image which can 
then be coded so that the different depth information is 
apparent in the image. 

In the preferred embodiment of the invention, the 
35 intensity variance data is processed in accordance with 
the phase image production method disclosed in our co- 
pending International Patent Application Nos . 
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PCT/AU99/00949 and PCT/AU02/01398 . The contents of these 
two International applications are incorporated into this 
specification by this reference. 

5 Thus, the step of processing the intensity variance data 
preferably comprises: 

(a) producing a representative measure of the 
rate of change of variance data values of said radiation 
wave field over a selected surface extending generally 

10 across the wave field; 

(b) producing a representative measure of said 
radiation wave field relating to the scene over said 
selected surface; 

(c) transforming said measure of rate of change 
15 of variance to produce a first integral transform 

representation and applying to said first integral 
transform representation a first filter corresponding to 
the inversion of a first differential operator reflected 
in said measure of rate of change of variance to produce a 
20 first modified integral transfoirm representation; 

(d) applying an inverse of said first integral 
transform to said first modified integral transform 
representation to produce an un transformed representation; 

(e) applying a correction based on said measure 
25 over said selected surface to said un transformed 

representation ; 

(f ) transforming the corrected un transformed 
representation to produce a second integral transform 
representation and applying to said second integral 

30 transform representation a second filter corresponding to 
the inversion of a second differential operator reflected 
in the corrected untransf ormed representation to produce a 
second modified integral transform representation; 

(9) applying an inverse of said second integral 

35 transform to said second modified integral transform 
representation to produce a measure of phase of said 
radiation wave field across said selected plane so as to 
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produce said image data as phase image data containing the 
depth information. 

The step of producing a representative measure of said 
5 radiation wavefield relating to the scene over the 

selected surface may use intensity values to obtain the 
representative measure or variance values . 

In particular, in one embodiment, intensity values at the 
10 selected surface are used. 

In another embodiment, values are taken representing 
maximum focus from any of the intensity data values. 

15 In another embodiment, variance data values are used as a 
representative measure . 

In a still further embodiment, maximum variance values 
taken from the sets of variance data values are used. 

20 

Preferably at least one of the first or second 
differential operator has a form based on an optical 
system used to acquire the radiation for producing the 
representative measure of the rate of change of variance 
25 of the radiation wavefield over the selected surface 
extending generally across the wavefield. 

Preferably both the first and second differential 
operators have a form based on the optical system. 

30 

PreferedDly the first and second integral transforms are 
produced using a Fourier transform. 

In one embodiment of the invention the differential 
35 operators have the form: 
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where , 

T, {p) = 2nidz jrjT^'^ (p, Tj)drj 
and 
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and wherein p is radial lateral spatial frequency, the r| 
longitudinal spatial frequency and & is the defocus 
distance in the plane of the object. Also 



NA 

^ ^ ^ ^ ^condensor 

~ NA 

^-'^ objective 

Where A^^,,„^,„,,, and A^^,,,,,,,,, are respectively the numerical 

aperture of the condenser and the objective (These are 
settings and dimensions on the microscope) . p^^^ is the 

20 maximum spatial frequency accepted by the objective. 



Preferably the step of taking at least two images of a 
scene comprises taking the first image at a first 
defocused plane to produce the first image as negatively 
25 focused image data and taking the second image of the 
scene at a second defocused plane to produce positive 
defocused image data, the negative and positive defocused 
image being taken on respective sides of a focal plane 
which would produce an in focus image of the scene. 

30 

Preferably the step of determining an intensity variance 
comprises determining a measure of the sharpness of each 
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of the dat:a values relative to the sharpness of data 
values surrounding that data value. 



10 



In the preferred embodiment of the invention the images 
are captured by a charge coupled device and an intensity 
value of pixels in the image is determined, and that 
intensity value is compared with the intensity of 
surrounding pixels in order to provide a variance value at 
each pixel . 

Preferably the variance is obtained by the following 
equation for each pixel: 



15 



wherein n is the particular pixel, and the values 1 to 9 
represent the pixels which surround that pixel in an array 
of pixels . 



Thus, in this embodiment, each pixel is compared with 
20 pixels which immediately surround the particular pixel in 
order to determine the variance . In other embodiments , 
each pixel could be compared with the pixels which 
surround that pixel, together with the pixels which 
surround the pixel and are one further pixel away and so 
25 on. The variance could also be determined by other 

methods so that a variance of each pixel relative to the 
other pixels in each data set is obtained. 



Preferably the step of coding parts of the image data 
30 comprises applying different colours to parts of the image 
which have the same depth information so as to produce a 
visual image in which the relative distance of parts of 
the scene compared to one another can be determined. 
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Preferably a third image of the scene is taken at an In 
focus plane to produce a third set of Intensity data 
values and the third set of Intensity data values are 
5 overlaid with the coded Image data to produce a visual 
Image containing structural Information of the scene as 
well as the different depth Information In the scene. 

However, In other embodiments, rather than providing a 
10 visual Image, the Image could be an electronic Image 

stored In a computer, and the code reference could simply 
be a code value allocated by the computer to the different 
depth Information so that If Images at a particular 
location are required. Information as to those parts of 
15 the Image can be extracted by the computer by the code 
reference . 

Preferably the different depth Information is provided by 
allocated a grey scale value to pixel values in the image. 

20 

The invention may also be said to reside in a computer 
program for performing the above-mentioned method. 

The invention may also be said to reside in an apparatus 
25 for producing an image containing depth information 
comprising : 

a Ccuaera for detecting a radiation wavef ield 
emanating from a scene and taking at least two images of 
the scene at different planes relative to the scene to 

30 produce a first set of intensity data values and a second 
set of intensity data values; 

a processor for determining an intensity 
variation of data values in the first set of data values 
compared to other data values in the first set of data 

35 values to produce a first set of variance data values, and 
for determining a variance of one data value in the second 
set of data values compared to other data values in the 
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second set of data values to produce a second set of 
variance data values ; 

the processor also being for processing the first 
and second variance data values to produce image data of 
5 the scene containing depth information; and 

the processor also being for coding parts of the 
image which have the same depth information with a code 
reference to identify the different depth information in 
the image . 

10 

Thus , the processor for determining an intensity variation 
is a processor for: 

(a) producing a representative measure of the 
rate of change of variance data values of said radiation 

15 wave field over a selected surface extending generally 
across the wave field; 

(b) producing a representative measure of said 
radiation wave field relating to the scene over said 
selected surface; 

20 (c) transforming said measure of rate of change 

of variance to produce a first integral transform 
representation and applying to said first integral 
transform representation a first filter corresponding to 
the inversion of a first differential operator reflected 

25 in said measure of rate of change of variance to produce a 
first modified integral transform representation; 

(d) applying an inverse of said first integral 
transform to said first modified integral transform 
representation to produce an un transformed representation; 

30 (e) applying a correction based on said measure 

over said selected surface to said untransf ormed 
representation ; 

(f ) transforming the corrected untransf ormed 
representation to produce a second integral transform 

35 representation and applying to said second integral 

transform representation a second filter corresponding to 
the inversion of a second differential operator reflected 
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in the corrected untransf ormed representation to produce a 
second modified integral transform representation; 

(g) applying an inverse of said second integral 
transform to said second modified integral transform 
5 representation to produce a measure of phase of said 

radiation wave field across said selected plane so as to 
produce said image data as phase dLmage data containing the 
depth information. 

10 The step of producing a representative measure of said 
radiation wavefield relating to the scene over the 
selected surface may use intensity values to obtain the 
representative measure or variance values . 

15 In particular, in one embodiment, intensity values at the 
selected surface are used. 

In another embodiment, values are taken representing 
maximum focus from any of the intensity data values. 

20 

In another embodiment, variance data values are used as a 
representative measure . 

In a still further embodiment, maximum variance values 
25 taken from the sets of variance data values are used. 

Preferably at least one of the first or second 
differential operator has a form based on an optical 
system used to acquire the radiation for producing the 
30 representative measure of the rate of change of variance 
of the radiation wavefield over the selected surface 
extending generally across the wavefield. 

Preferably both the first and second differential 
35 operators have a form based on the optical system. 

Preferably the first and second integral transforms are 
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produced using a Fourier transform. 

In one embodiment of the invention the differential 
operators have the form: 

Tp 

where , 
and 



15 



and wherein p is radial lateral spa-bial frequency, 'the x\ 
longitudinal spatial frequency and dz is -the defocus 
distance in the plane of the object. Also 



20 



condemor 



NA 



objective 



25 



where NA^^^^^^^^ and NA^^j^^^^^^ are respectively the numerical 

aperture of the condenser and the objective (These are 
settings and dimensions on a microscope, if a microscope 
is used in the capture of the various images by camera) . 
p^fj^j is the maximum spatial frequency accepted by the 

objective . 



If images are captured through a different optical system 
30 such as a telescope or simply by the optics of a camera, 
then the above operators are modified to suit the various 
optics of that system. 
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10 



15 



Preferably the step of taking at least two images of a 
scene comprises taking a first image at a first defocused 
plane to produce a negatively focused image and taking a 
second image of the scene at a second defocused plane to 
produce a positive defocused image, the negative and 
positive defocused image being taken on respective sides 
of a focal plane which would produce an in focus image of 
the scene . 

Preferably the step of determining an intensity variance 
comprises determining a measure of the sharpness of each 
of the data values relative to the sharpness of data 
values surrounding that data value. 



In the preferred embodiment of the invention the images 
are captured by a charge coupled device and an intensity 
value of pixels in the image is determined, and that 
intensity value is compared with the intensity of 
20 surrounding pixels in order to provide a variance at each 
pixel . 

Preferably the variance is obtained by the following 
equation for each pixel: 

25 



wherein n is the particular pixel, and the values 1 to 9 
represent the pixels which surround that pixel in an array 
of pixels and I is the intensity value of the respective 
30 pixel. 



PrefercJ:>ly the step of coding parts of the processed image 
comprises applying different colours to parts of the image 
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whxch have the same depth information so as to produce a 
visual image in which the relative distance o£ parts of 
the scene compared to one another can be determined. 

5 However, in other embodiments, rather than providing a 
visual image, the image could be an electronic image 
stored in a computer, and the code reference could simply 
be a code value allocated by the computer to the different 
depth information so that if images at a particular 
10 location are required, information as to those parts of 

the image can be extracted by the computer by reference to 
the code reference. 

Preferably the different depth information is provided by 
15 allocated a grey scale value to pixel values in the image. 

Brief Description of the Drawings 

A preferred embodiment of the invention will be described 
with reference to the accompanying drawings in which: 
20 Figure 1 is a diagram to assist in illustrating 

the preferred embodiment of the invention; 

Figure 2 is a view of an apparatus embodying the 

invention ; 

Figure 3 is a diagram showing pixels of a Ccuaera 
25 used to take the images in the preferred embodiment of the 
invention ; 

Figure 4 is a diagram illustrating the production 
of different depth information; 

Figure 5 is a diagram illustrating how different 
30 depth information is coded according to the preferred 
embodiment of the invention; 

Figure 6 is a flowchart showing an implementation 
of the method of the preferred embodiment; 

Figure 7 is a schematic illustration of an 
35 arrangement which uses a microscope in the capture of 
images ; 

Figure 8 is a schematic drawing of an exemplary 
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system according to one embodiment of the invention; and 

Figure 9 is a flow chart relating to one 
embodiment of the invention. 



5 Detailed Description of the Preferred Embodiment 

The preferred embodiment of the invention effectively 
processes images in accordance with the phase 
determination technique disclosed in our International 
Applications PCT/AU99/00949 and PCT/AU02/01398 . However, 

10 before the algorithm disclosed in these applications is 
used to process the data relating to the captured scenes 
to produce the image, the data is first manipulated so as 
to produce a variance data so the variance data is 
processed in accordance with the algorithms rather than 

15 pure intensity data, as disclosed in the above 

International applications. However, apart from this, the 
method of processing the data is the same as that 
disclosed in the above International applications. 

2 0 Figure 1 is a diagram showing an image of a scene at three 
different focal planes. The plane Iq is an in focus plane 
which will be used to capture an in focus image of the 
scene . The scene may include a niimber of articles or 
objects such as vehicles, personnel, or the like, and in 

25 general, the preferred embodiment of the invention is to 
be used for obtaining thermal images . However , the 
invention could be used with other types of images and 
improve the depth information in those dLmages . The 
invention has particular application to thermal images 

30 because thermal images generally contain little or no 
depth information which enables different parts of the 
images to be compared in teirms of distance with other 
parts of the image. Thus, in thermal images, it is very 
difficult to ascertain whether one part of the image is 

35 closer to or further away from the location from where the 
images were captured, or where parts of the image are 
located in space relative to one another. 
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The plane V- is a negatively defocused plane within the 
depth of field of the camera 20 (see Figure 2) used to 
capture the images, the plane Iq is the focal plane, and 
5 the plane V+ is a positively defocused plane also within 
the depth of field of the camera 20, but positively 
defocused relative to the plane Iq. In other words, the 
plane V- is closer to the scene than the in focus plane Iq 
and the plane V+ is further away from the scene than the in 
10 focus plane Iq. 

Typically, parts of the image taken at the plane V. which 
happen to be in focus t the plane V- will produce a very 
sharp image point 21 at the plane but will be a very 

15 fuzzy defocused image at the plane V+ as shown by reference 
22 . The waveform which produces the images at 21 and 22 
is shown at 23 in Figure 1. Similarly, part of the image 
24 which is defocused and blurry at the plane V- may be 
focused at the plane V+ as represented by reference 25 . 

2 0 Once again, the form of the wavefield producing these 
parts of the images is shown at 26. 

As is shown in Figure 2, the camera 20 is used to capture 
images of the scene schematically shown at S in Figure 2 
25 at the planes Iq and V+ so that three sets of intensity 

data at those three planes is captured by the pixels of 
the camera 20. 

In the preferred embodiment of the invention, the camera 
30 20 is a cryogenically-cooled, 640 x 480 pixel array 

thermal camera operating in the 8 to 12 micron band with 
adjustable-focus 6e lens . The intensity data captured by 
the camera 20 is processed in a processor 40 and images 
may be displayed on a display 42 . 

35 

The processor 40 first determines a variance value for 
each of the pixels for each of the two images at the 
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planes V. and V+. Thus, two sets of variance data xs 
therefore produced. 

Figure 3 shows an example of nine pixels of the pixel 
5 array forming camera 20 and the variance at the pixel X5 Is 
determined In accordance with the following equation: 



10 Thus, a variance value for the fifth pixel In the array of 
Figure 3 Is determined. In Figure 3, pixels which are 
arranged at the edge of the array can be set to a value 0 
so that they are not Included In the calculations . This 
reduces the size of the Image slightly by 2 pixel widths. 

15 However, It makes calculation much easier because all of 
the variances can be detezrmlned In the above mentioned 
manner because all pixels will then be surrounded by 8 
other pixels. Similarly, variance values for each of the 
other pixels in the array is determined for each of the 

2 0 image pixels of planes V. and V+ so that two sets of 

variance data is produced. The two sets of variance data 
is then used in the algorithm to produce the phase image 
data of the scene instead of the pure intensity data 
originally captured by the camera 20. The use of the 

25 variance data provides a measure of the sharpness of the 
image at each of the planes V- and V+ and therefore, 
inherently Includes a measure of the depth of parts of the 
Image relative to other parts of the image. The variance 
data is then processed in the same manner as the intensity 

30 data, as disclosed in the aforementioned International 

applications , so as to produce a phase Image of the scene 
in which depth of field information Is Included. Because 
the variance data Is used to produce the phase image, the 



1 

' 8 
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image will contain very little detail of the actual 
structure of articles or things in the scene, but will 
contain general outlines showing features of the scene. 
The depth of position information in the processed image 
5 can be given by a grey scale value at each of the pixels 
in the processed image data, which in turn therefore 
provide a relative measure of the part of the scene which 
relates to each pixel relative to other parts of the 
scene, as well as relative to the position of parts of the 
10 scene from where the images were taken. For example, as 
is best shown in Figure 4, something which is very close 
to where the image is taken appears darker, as shown by 
reference 61, and something which is further away appears 
much lighter as at 62 . 

15 

The grey scale values allocated to each of the pixels in 
the processed phase image can be used to provide an actual 
measure of the distance of the parts of the scene from 
where the image was taken and relative to one another. 

2 0 For example, if the optics of the camera is known and the 
distance between the planes Vi and V+ is also known, actual 
measurement values can be associated with various grey 
scale values in the image to provide a distance measure of 
parts of the scene relative to one another, and also from 

25 where the images were taken. 

The image shown in Figure 4 can then be coded so as to 
more clearly bring out the different distances of articles 
in the image by applying a code reference to parts of the 

30 image which have the same grey scale value. This can be 
done in a variety of different ways, but most preferedDly 
is performed by coding parts of the image with a false 
colour palette so that parts of the image which have the 
same grey scale value are coloured with the same colour. 

35 This therefore makes it very much easier to perceive depth 
information in the image than in the mere conventional 
thermal image in which the entire image appears to be 
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approximately the same shade. For instance, in Figure 5, 
the grey scale values are shown for near articles in the 
image which may be say 10 metres from where the image is 
taken, up to 100 metres away- Those grey scale values can 
5 be allocated with different colours in the manner referred 
to above, so that a particular colour can be associated 
with a particular distance in the image to make it 
therefore much more easy to perceive where articles in the 
image are actually located relative to one another. 

10 

In order to produce an image in which the details of the 
article are shown, rather than the pure distance 
information, the processed image can then be applied to 
the intensity image taken at plane lo which is the in focus 

15 plane, to give both the intensity information of the image 
and the depth information of the image in the one image. 
Thus, the intensity information provides the structural 
detail of the various components in the image and the 
depth information, as previously described, enables you to 

2 0 ascertain the relative position of parts of the image. 
Thus, this image will not only contain the depth 
information, but will also make it much easier to 
ascertain what various components of the image are, 
because the actual structure in the image will be seen. 

25 

Thus, according to the preferred embodiment of the 
invention, the imaging range of the lens on the camera 20 
and the colour rendering can be calibrated and therefore 
positions of the image can be quantified. This technique 

30 then offers the opportunity to provide a passive ranging 
capability. In conditions of near thermal uniformity, 
where there is very little contrast in the field of view, 
the depth positional information provided by this process 
suggests another contrast mechanism for viewing and 

35 processing images of scenes. Finally, under obscured 

viewing conditions , such as under camouflage , the outline 
of an object can be disrupted, making identification 
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difficult:. The colour depth rendering of the scene allows 
the piecing together of image components at common 
positions, as even though the outlines can be disguised, 
the physical positions in space cannot be. 

5 

The manner in which the algorithm operates to process the 
variance data to produce the phase image is described 
below . 

10 At each point in space, an optical beam possesses two 

properties: Intensity and phase. Intensity is a measure 
of the amount of energy flowing through each point, while 
phase gives a measure of the direction of the energy flow. 

15 Intensity may be measured directly, for example by 

recording an image on film. Phase is typically measured 
using interference with a "reference beam" . In contrast 
the present method gives a non-interf erometric method for 
measuring phase. 

20 

Intensity can be measured over two parallel planes V., V+ 
extending across the direction of propagation of the wave 
field on the side remote from the incident radiation. 

25 The present invention determines phase by providing a 
solution to the transport-of -intensity equation: 

(1) V^.(/V^^ )=-*^ 

30 where I is the intensity in the plane, the gradient 

operator in the plane is denoted , k is the wave number 

of the radiation, and dl/dz is the intensity derivative or 
rate of change of intensity. However, in order to obtain 
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"the dep'th Information which is previously referred to, 
rather than using the intensity values which are obtained 
by the camera at the two planes , the variance values 
calculated in the above mentioned manner are obtained so 
5 that a variance derivative, or rate of change of variance 
dVldz . Note that dVldz is estimated from the difference of 
the measurements in the planes V. and V+ shown in Figure 1, 
while the intensity I is given by the average of the 
measurements . 

10 

In order to obtain a solution to equation 1 the function A 
is first defined as: 

(2) Vj,^ = FV^^. 

15 Thus equation (1) becomes: 

(3) V^#(V^^) = -^a,F . 

Making use of a standard identity Vj^ • Vj^ = Vj^^ , this may be 
2 0 written : 

(4) V^A^-kd^V 

where Vj^^ denotes the two-dimensional Laplacian acting 
2 5 over the surface of the image. This equation has the 
following symbolic solution: 

(5) ^ = -A:V/^5,F . 

30 If the gradient operator Vji^ is applied to both sides of 
this equation, it becomes: 
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(6) V^^ = -A:V^V/Xf . 

The defxning equation (2) for the function A allows (6) to 
5 be transformed into: 

(7) rs/j = -kv^'^^-''d^v . 

Dividing both sides by I then yields : 

10 

(8) V_,^ = -^rV_LV/'5,F . 

However, as is mentioned above, rather than dividing both 
sides by the intensity value I, a maximum focus value 
15 across all of the images could be used, or the variance 
data at the plane could be used, or a maximum variance 
across all of the variance state of values could be used. 

Taking the two dimensional divergence Vj^ • of both sides of 
20 (8) , and again making use of the standard identity 

Vj^ • Vj^ = Vj^^ , then (8) becomes : 

(9) V = -k^^ • [/'V^V^-'5,fJ . 

2 5 This equation has the following symbolic solution: 

(10) = -wi(v^ • [rV^v/'a.rJ) . 

In order to implement a practical solution to equation 
30 (10) , the following formulae are required. A suitaJDly- 
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well-behaved function f (x,y) may be written in the form of 
a two-dimensional Fourier integral: 

(11) f{x,y)= J J hK^ky'''''''dk^dk^. 

5 

The function f{k^,ky) is called the ^^Fourier transform" of 
f (x,y) . 



10 



The X derivative of (11) yields; 



(12) -^Ax^y) = J J [ ik^ K K,K) dk^dk^. 



Hence the Fourier transform of — f (x,y) is equal to the 

ax 

Fourier transform of f(x,y) multiplied by ik^ , Stated 

15 differently, — = zF'^Xr^F , where F denotes Fourier 

dx 

transf oinnation and F~^ denotes inverse Fourier 
transformation. Similar considerations apply to 

(x,y) . 

dy 

If the Laplacian = — t'^^TT obtained and 

20 similar reasoning applied, it follows that V^^ = -F'^fc^ , 
where =^kl +kl . Thus : 
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Here, F denotes Fourier transformation, F"^ denotes inverse 
Fourier transformation, (k, ky) are the Fourier variables 
conjugate to (x,y) , and 



Equations (13) can be used to rewrite equation (10) in the 
form 



10 



(14) 



dV 
dz 



dz 



In practice division by intensity is only performed if 
that intensity is greater than a certain threshold value 
15 (eg. 0«1% of the maximum value) . 

Division by kr does not take place at the point kr = 0 of 
Fourier space; 

instead multiplication by zero takes place at this point. 
20 This amounts to taking the Cauchy principal value of the 

—2 

integral operator Vj^ 



In order to quantitatively measure the phase of object it 
is necessary to incorporate some physical constants into 
25 the phase recovery algorithm given in Equation (14) 

relating to the experimental setup in use to quantify the 
variables kx, ky. This can be done by rewriting equation 
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(14) In -the following form suitable for ±mplemeni:a'b±on 
using a fast Fourier transform: 



XSz{N^xf /(x,y) /^+/^^ 

XSz(NAxf i^+f /(x,y) e+f 



where /^ye 



-TV 



Index the frequent components of 



2 2 

where the Intensity derivative dy{x^y)Ls obtained by 
subtracting two Images and VL separated by a distance 
Sz , i and j are the pixel numbers on the Image ^ and using 
10 the fact that the Fourier space step size Is given by 



where the Image Is an N x N array of pixels of size Ax . 

15 Thus In addition to measuring the three Intensity 

distributions It Is necessary to know the pixel size Ax , 
defocus distance dz and wavelength X In order to make a 
quantitative phase measurement. All of these quantities 
can be readily determined: the pixel size can be 

2 0 determined directly for example from the CCD detector 

geometry (In the case of direct Imaging) , or by existing 
techniques for calibrating the transverse distance scales 
(In the case of an Imaging system) , the defocus distance 
can be measured directly, and the spectral distribution of 

2 5 the Illumination can be determined either by 

monochromatlng the Incident field or by analysing the 
spectral distribution of the radiation using existing 
spectroscopic methods . 
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An example of the phase-retrieval method Implementing the 
solution of equation (14) can be represented by the 
flowchart shown in Figure 6. As shown in Figure 2 the 
5 quantitative determination of phase of a radiation wave 
field proceeds from a set of variance measurements {V„} 
over the two spaced apart planes V- and V+. A measurement 
of central intensity I(x,y) in a selected plane parallel 
to and midway between the planes V- and V+ is also 

10 obtained. The intensity measurements are performed over a 
defined array on each of the two planes A and B and the 
respective values subtracted to produce a measure of the 
intensity derivative. This value is multiplied by the 
negative of the average wave number. The data are split 

15 into two component sets and a fast Fourier transform is 

performed to produce the respective x and y components in 
the Fourier domain. A filter is then applied to the 
Fourier doiaain representations to correspond to the 
inversion of a differential operator reflected in the 

20 un transformed representation. An inverse Fourier 
transform is then performed on each of the x and y 
coxnponents to produce a spatial domain value from which 
the differential operator has been removed. A division by 
the central intensity I(x,y) obtained by averaging the 

25 intensity measurements over planes V+ and V- is then 
performed if the intensity level is above a selected 
threshold level . The resultant data is again Fourier 
transformed and multiplied by the same filter to again 
correspond to the inversion of a differential operator 

30 reflected in the un transformed data. The resultant 

components are again inverse Fourier transformed and 
summed to provide a retrieved phase measurement. 
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will be apparent that in general the method according 
to this invention can proceed from any suitable 
representative determination of intensity derivative or 
rate of change of intensity over a selected surface 
5 extending across the propagation direction and the 

intensity over that same surface. As will be explained in 
various exeunples these data can be obtained in a variety 
of ways and the method implemented to yield phase of the 
radiation wave field. 



10 



Rewriting equation (14) with: 



-2 



15 



0,y{k^,ky,a) 



kyk^ 



-2 



gives 



20 



(15) 




where : 



25 



^(x^j;) denotes the recovered phase ^ 



F denotes Fourier transformation, and F 



denotes 



inverse Fourier transformation. 
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^(x,y) is the Intensity distribution over the plane 
of Interest, 

(x,y) are Cartesian coordinates over the plane of 
interest, 

5 (k^.ky) are the Fourier variables conjugate to (x,y) 

k=27t/x is the average wave number of the 
radiation , 

/I is the average wavelength of the radiation, 
dV/dz is the estimate for the longitudinal variance 
1 0 derivative , 

a is the regularization parameter used to 
stabilize the algorithm when noise is present. 

As given above, the solution to the transport of intensity 
15 equation (1) assumes a perfect imaging system. That is, 
there are no ^^aberrations" present in the optical system 
used to obtain the intensity data which is fed into the 
algorithm. Of course, no imaging system is perfect. The 
imperfections present in an imaging system may be 
2 0 quantified by a set of numbers : 

(16) ^,,^2.^3 — 

which are termed aberration coefficients . 

25 

If intensity data were taken on an imperfect instrximent 
whose imperfections were characterized by a certain set of 
known aberration coefficients ^,,^29^3'—/ would be 
desirable if the filters Cl^(k^ky^a) and Q.y{k^ky^a) present in 

30 (15) could be replaced by modified filters which 

explicitly depend upon the aberration coefficients : 
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(17) n^(k^ky,a,A^,A2,A^,.,.) and Cly{k^ky,a,AiM2M3^") 

Thxs would allow the xmperfections o£ the xmaglng system 
5 to be explicitly taken into account, leading to 

quantitatively correct phase retrieval using aberrated 
imaging systems . For the special case of a non-absorbing 
phase object in a radiation wave field of uniform 
intensity with weak (i.e. much less than 2;r radians) phase 
10 variations the appropriate modified filters lead to the 
following functional form for the phase-retrieval 
algorithm: 



(18) iP(x,y) = F' 

1 5 where : 



20 



25 



_F {Kterrated(X,y)-l} 



Vkberrated (x ^y) is the aberrated variance measured at 
def ocus distance dz , 

Am are the aberration coefficients which 
characterize the ixiKperfect imaging system. 

If a filter is defined: 



(19) 0.{k,ky,aMx,A2,A^,...) 

1 

4- 2n,dzMkl +kl)-2Z^ A^^k^k; 



Then (18) becomes : 

(20) Hx>y) = F'nFlF'nF{\,^^(x,y)-l } 
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The term { I aberrated (^^y)~^ }is a measure of rate of change of 
intensity. Iq Intensity is a measurable constant for 
uniform intensity so that (20) is the Scune general form as 
(15) . Consequently the special case of aberration can be 
5 dealt with by changing the filter in the general method 
described above. The x and y component filters and 
are given by 

(21) 

10 

Figure 7 schematically shows an arrangement for 
q[uantitative phase amplitude microscopy. A sample is 
illuminated using a source of white light Kdhler 
illumination 15, commonly found on optical microscopes. 

15 The light is transmitted through an object 16 and 

collected by a microscope imaging system 17 and relayed to 
a CCD camera 18 or other digital imaging device having a 
planar imaging surface. Three images are collected: an 
in-focus image, lo, and two slightly out of focus images I^. 

20 and The defocus is obtained by suitable means such as a 
drive system 19 to adjust the microscope focus knob. The 
defocus introduced is usually quite small so that 
degradation in spatial resolution is minimised, although 
the optimal amount of defocus to use is determined by 

25 sample properties and imaging geometry such as 
magnification , numerical apertures , etc . 

When taking the images the niamerical aperture of the 
condenser is chosen to be less than the numerical aperture 
30 of the objective being used. If this is not the case then 
serious image degradation will occur, although the precise 
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amount by which t:he condenser and objective numerical 
apertures should differ involves a tradeoff between image 
fidelity and spatial resolution, with the optimal 
difference depending on the Scuaple properties and the 
5 optics used. 

The variance data values are determined in the manner 
previously described from the intensity values collected 
at the planes V+ and V. and are subtracted to produce a 
10 representative measure of rate of change of variance 

(variance derivative) • From this and the intensity data 
of collected image lo the method described above can be 
used to produce quantitative information about the 
magnitude of the phase shift in the image plane. 

15 

There may be cases in which it is desirable to take more 
than two images in order to obtain a better estimate of 
the intensity derivative dV/dz . A function can then be 
fitted to this data from which dV/dz can be computed and 
2 0 used in the phase determination method in place of the 
simple subtraction of two images normally used. 

It is also possible to operate this system in reflection 
geometry to obtain surface topography. The principle of 
25 operation is the same, however the optics have to be 

folded back on themselves to form a reflection geometry - 
otherwise the process is identical . 

For certain applications it can also be desirable to 
30 filter the light to a particular wavelength, although this 
is not necessary for the described imaging process as it 
works equally well with white light. 
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An implementation is shown in Figure 8. An Olympus BX-60 
optical microscope 20 was equipped with a set of UMPlan 
metallurgical objectives and a universal condenser to 
provide Kohler illumination. In order to be able to 
5 compare the results with existing imaging modes Nomarski 
Die optics and a set of cover-slip corrected UplanApo 
objectives were also acquired for this microscope, 
enabling images to be taken of the same field of view 
using both phase retrieval and Nomarski DIG for the 
10 purposes of qualitative comparison. A 12 -bit scientific 
grade Photometries SenSys CCD Ccuaera 21 equipped with a 
1300x1035 pixel Kodak KAF-1400 CCD chip was added to the 
0 . 5x video port on the microscope to acquire digital 
images of the sample. 

15 

The phase recovery technique of this embodiment of the 
invention requires the acquisition of def ocused images . A 
stepper motor drive system 22 was attached to the focus 
knob of the microscope. This stepper motor 22 was 

20 attached to the parallel port of a 133MHz Pentium PC 23 
also used to control the CCD camera 21, enabling full 
automation of the acquisition of through-focus image 
sequences. This data acquisition system was linked to 
custom software written to recover phase images from the 

2 5 CCD images, thereby enabling full automation of the image 
acquisition and data processing sequences . 

The form of the differential operators used in the 
preferred embodiment of this invention are based on the 
30 optics of the system used to obtain the above-mentioned 

images . Thus , the operator takes into account the details 
of the optical system used to take the images. This is 
achieved by: 
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Det.ermi.ne 'the numerical aperture o£ the condenser^ NAcondensor 

Determine NA of objective, NAobjectivc , and pobjectivef the radius 
5 of the objective aperture 

NA 

~ NA 

-^^ "^objective 

(These are settings and dimensions on the microscope . ) 

10 

Determine radial lateral spatial frequency, p, and 
longitudinal spatial frequency, T| . 

(These are dependent on the pixelation and position 
15 distribution of images taken in the series.) 

Determine the wavelength, X, of the radiation to be used. 
The new form of the operator is 

20 




where , 

25 7> (p) = 2m Sz jrjT^'^ (p, Tj)drj 
and 




30 

Figure 9 is a flow chart generally Illustrating how Tp is 
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determined by means of the above eG[uatlon merely showing 
breakdown of the various components of the equation. 

The calculator depth image can then be presented in 
5 different modalities to give better visualisation. Such 
modalities including DIG, bright field images, dark field 
images and other conventional modality images. Techniques 
for forming these types of images from the phase data 
determined in accordance with the present invention are 
10 described in our copending International Application No. 

PCT/AU02/00590 (the contents of which is incorporated into 
this specification by this reference) . 

Since modifications within the spirit and scope of the 
15 invention may readily be effected by persons skilled 

within the art, it is to be understood that this invention 
is not limited to the particular embodiments described by 
way of exaunple hereinabove. 

20 In the claims which follow and in the preceding 

description of the invention, except where the context 
requires otherwise due to express language or necessary 
implication, the word ^^comprise'' , or variations such as 
^^comprises" or ^^comprising" , is used in an inclusive 

2 5 sense, ie. to specify the presence of the stated features 
but not to preclude the presence or addition of further 
features in various embodiments of the invention. 



